Stabilization of finite element method
When performing calculations using the finite element method, we perform following steps
- We determine the area where we want to perform our calculations and generate a grid of elements covering a given area. For example we consider the area \( \Omega=[0,1]^2 \) covered with a mesh of 4 by 4 elements (patch of 16 elements on which we will later span the B-spline base functions).
- We choose a differential equation in the strong form that models the physical phenomenon that we want to simulate in our area, for example the heat-transfer problem \( -\Delta u=f \).
- We establish the boundary conditions that describe what happens to the boundary of our area, for example \( u=0 \textrm{ for } x=0,y\in(0,1); u=1 \textrm{ for } x=1,y\in(0,1), \frac{\partial u}{\partial y}=0 \textrm{ for }x\in(0,1),y\in \{ 0,1 \} \).
- We transform the strong formulation into the so-called weak formulation, for example \( b(u,v)=l(v); b(u,v)=\int_{\Omega} \nabla u \cdot \nabla v dxdy; l(v) = \int_{\Omega} x v(x,y) dxdy \) where the x function on the right-hand side is the extension of the boundary condition into the entire area (i.e., such a function that takes 0 on the left boundary of the area and 1 on the right and her derivatives in the y direction are equal to 0 as required by the boundary conditions).
- We will establish a finite-dimensional base of functions used for approximation of our solution, for example using the knot vectors on ours element grid [0 0 0 1 2 3 4 4 4] on the x axis and [0 0 0 1 2 3 4 4 4 ] on the y axis, which corresponds to the basis of the B-spline functions \( B^x_{i,2}B^y_{j,2}, i,j=1,...,6 \). They create approximation space \( U_h = span \{ B^x_{i,2}B^y_{j,2} \}_{i,j=1,...,6} \).
- We will establish a finite-dimensional basis of functions for testing of the strong formulation, for example the same as for approximation, using knot vectors on our mesh of elements [0 0 0 1 2 3 4 4 4 ] on the x-axis and [0 0 0 1 2 3 4 4 4] on the y-axis, which corresponds to the base B-spline functions \( B^x_{k,2}B^y_{l,2}k,l=1,...,6 \). They form the testing space \( V_h = span \{ B^x_{k,2}B^y_{l,2} \}_{k,l=1,...,6} \).
- Thus, we create a discrete version of our problem \( b_h(u_h,v_h)=l(v_h) \forall v_h \in V_h \), where in we replace the abstract perfect solution \( u \) by the approximate solution \( u_h = \sum_{i,j=1,...,6} u_{i,j} B^x_{i,2}(x)B^y_{j,2}(x) \), and as testing functions we choose the elements of a finite-dimensional testing space \( v_h = \sum_{k,l=1,...,6} v_{k,l} B^x_{k,2}(x)B^y_{l,2}(x) \).
- Generate a system of linear equations by inserting into the weak formulation our base functions for approximation and for testing \( \{ \int_{\Omega} \frac{\partial B^x_{i,2} } {\partial x } B^y_{j,2} \frac{\partial B^x_{k,2} }{\partial x} B^y_{l,2}+ B^x_{i,2} \frac{\partial B^y_{j,2} }{\partial y } B^x_{k,2} \frac{\partial B^y_{l,2} }{\partial y } \}_{i,j=1,...,6;k,l=1,...,6 } \{ u_{k,l} \}_{k,l=1,...,6} = \{ \int_{\Omega} x B^x_{k,2}B^y_{l,2} \}_{k,l=1,...,6} \)
- We solve the obtained system of linear equations, using for example LU factorization algorithm with pivoting.
- We modify our solution by subtracting the shifted condition \( u_h = u_h-x \).
This chapter focuses on the issue of stabilization of the finite element computations. The most important mathematical theorem allowing the study of the stability of the finite element method was independently proposed by prof. Ivo Babuśka, prof. Franco Brezzi, and prof. Olga Ładyżenska [1], [2].
They discovered (at about the same time independently of each other) equivalent conditions, now called inf-sup ("inf-sup condition"). The inf-sup condition is still used to this day by scientists who study the stability of the finite element method. This condition can be defined in abstract infinite dimensional mathematical spaces \( \inf_{u \in U } \sup_{v \in V } \frac{b_h(u,v)}{\|u\| \|v\| } = \gamma > 0 \) or in finite dimensional spaces resulting from our approximation of the solution with the basis of the B-spline function, used to approximate and test the solution \( \inf_{u \in U } \sup_{v \in V } \frac{b_h(u,v)}{\|u\| \|v\| } = \gamma > 0 \). It may happen that even if mathematicians have proved the fulfillment of the inf-sup condition over infinite dimensional spaces - that is, when we will start to solve a problem in finitely spaces dimension - the inf-sup condition will no longer be met. \( \inf_{u \in U } \sup_{v_h \in V_h } \frac{b_h(u,v_h)}{\|u\| \|v_h\| } = \gamma > 0 \) then even if the approximation space were infinitely dimensional, our simulation will be unstable anyway.
The problem with the inf-sup condition is that by selecting our finite dimensional sets (so-called bases) of functions, we do not have an approximation of the solution to test our equations to be sure that the solution obtained with their help on the computer will be sufficiently accurate. The supremum may no longer be realized when we restrict our test space to the finite dimensional space. To check this, we apply very much advanced math tools such as the inf-sup condition. For a large class of problems solved on computers using the finite elements method, this condition is satisfied, and we get quite exact solutions. However, there is also a fairly large class of problems for which the inf-sup condition is not met. This condition can be written as follows. We choose the version formulated by prof. Ivo Babuśka (all versions of this condition are described in the work of Leszek Demkowicz, ICES Report 0608, 2006 "Babuśka <=> Brezzi?" ) [3]:
\( \inf_{u_h \in U_h } \sup_{v_h \in V_h } \frac{b_h(u_h,v_h)}{\|u_h\| \|v_h\| } = \gamma_h > 0 \) If this condition is met, it will be possible to solve the system of equations generated using the finite element method and solving it will be a pretty good approximation to the ideal solution \( u \approx u_h = \sum_{i,j=1,...,6} u_{i,j} B^x_{i,2}(x)B^y_{j,2 }(x) \).
If this condition is not met, then either the generated layout equations will not be solvable for example, the Gaussian elimination algorithm will tip over because of the diagonal of the matrix will result in a numeric zero (a very small number which will not divide a row in the matrix), or the obtained solution of the system equations will be an incorrect solution to our problem, e.g. numerical oscillations will arise. In such a situation, it is necessary to stabilize our problem, which relies on:
- functional modification \( b_h(u_h,v_h) \) (on adding any additional stabilizing members to it, which is the case approximation and testing space of infinite dimensional, v ideally. These are zero, but in terms of space, they are finite dimensions which are not zero and improve the stability of our problem - an example of this method is the stabilization of the SUPG for equations advection-diffusion), or
- modifying the equations of the problem and the standard in which we try to meet inf-sup condition UhUh (for example in the Discontinuous Galerkin DG method, we break functions and modify the equations by adding terms equal to zero in infinite dimensional space, but not equal to zero in the approximation space. We modify the equation on the discrete level and the norm which gives the inf-sup condition also satisfied at the level of finite-dimensional discrete spaces), or
- on such reformulation of the problem that the approximation space and the test space are different and that the size of the space could be increased (towards the infinite dimensional test space where the supremum is met) \( V_h \) (on taking another space testing than of the approximation space, thanks to which the inf-sup condition will be a better approximation - this is how, for example, the method of minimization of the reziduum works).
All these stabilization methods make it possible to mathematically prove the inf-sup condition has been met. However, the mathematical tools needed to check that this condition is met are very mathematically advanced, so we will limit ourselves only to the application of several stabilization methods (SUPG, residuum minimization and the DG method) which will improve the stability, but we will not provide explanations as to why they work.
The method of minimizing the residuum also lies at the heart of the Discontinuous Petrov-Galerkin (DPG) method. This DPG method proposed by prof. Leszek Demkowicz is very popular in stabilization of adaptive finite elements computations. In the DPG method, the approximation and test spaces are broken to enable local static condensation on individual finite elements [4], [5], [6].
Bibliography
1. Ivo Babuśka: Error-Bounds for Finite Element Method, Numerische Mathematik 1971, dostęp:18.10.20192. Franco Brezzi: On the Existence, Uniqueness and Approximation of Saddle-Point Problems Arising form Lagrange Multipliers, R.A.I.R.O., Italy 1974, dostęp:18.10.2019
3. Leszek Demkowicz: Babuśka <=> Brezzi, The University of Texas at Austin, USA 2006, dostęp:18.10.2019
4. T.J.R. Hughes, L.P. Franca, M. Mallet: A new finite element formulation for fluid dynamics: VI. Convergence analysis of the generalized SUPG formulation for linear time-dependent multidimensional advective–diffusive systems, Computer Methods in Applied Mechanics and Engineering, Elsevier 1987, dostęp:18.10.2019
5. Daniele Antoni Di Pietro, Alexandre Ern: Mathematical Aspects of Discontinuous Galerkin Methods, Springer, Paris 2011, dostęp:18.10.2019
6. Leszek Demkowicz, Jay Gopalakrishnan: Discontinuous Petrov–Galerkin (DPG) Method. Part 2. Fundamentals, Encyclopedia of Computational Mechanics, Wiley 2017, dostęp:18.10.2019